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ABSTRACT 


A  Monte  Carlo  procedure  was  developed  to  evaluate  the  energy, 
angular  distribution,  and  intensity  of  either  the  scattered  neutron 
or  garama-ray  flux  that  penetrates  a  multibend  duct.  The  procedure 
has  been  coded  for  the  IBM  7090. 

A  detailed  presentation  of  the  Monte  Carlo  method  as  an 
approximation  to  the  Neumann  series  solution  of  the  Integral  trans¬ 
port  equation  Is  given.  Sampling  techniques  utilized  by  the 
procedure  are  described.  Included  in  these  techniques  are  splitting, 
Russian  Roulette,  statistical  estimation,  and  a  method  of  biasing 
the  sampling  from  the  source  angular  distributions. 

Results  obtained  with  the  procedure  are  compared  with  the 
data  taken  In  the  duct  penetration  and  systemlzatlon  experiment 
conducted  at  General  Eynamlcs/Port  Worth.  This  comparison  confirms 
the  validity  of  the  methods.  Further,  It  shows  that  the  procedure 
will  be  a  valuable  aid  In  the  analysis  of  experimental  data  as  well 
as  In  the  determination  of  the  validity  and  range  of  applicability 
of  some  of  the  simpler  methods  developed  to  calculate  the  flux 
penetrating  a  multibend  duct. 
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REPORT  SUMMARY 


Due  to  the  success  of  a  recent  application  of  the  Monte  Carlo 
method  to  problems  Involving  neutron  scattering  and  penetration 
through  straight  cylindrical  ducts,  it  was  decided  to  modify  the 
procedure  so  that  multibend  ducts  might  be  considered.  Routines 
providing  for  the  treatment  of  gamma-ray  scattering  and  for  the  use 
of  splitting  and  Russian  Roulette  techniques  were  also  added  to  the 
procedure.  One  of  the  greatest  assets  of  the  procedure  is  that  a 
rigorous  treatment  of  multiple  scattering  within  a  multibend-duct 
configuration  is  possible.  A  neutron  or  gamma-ray  source  may  be 
described  with  a  set  of  from  1  to  30  point  sources.  The  energy 
and  angular  distribution  and  the  intensity  of  the  scattered  flux 
are  calculated  for  each  of  a  set  of  from  1  to  30  detector  points. 
The  unscattered  flux  is  also  calculated  for  each  detector  position 
and  recorded  separately  so  that  the  total  flux  as  well  as  the  rela¬ 
tive  contributions  from  scattered  and  unscattered  radiation  may  be 
determined. 

The  Monte  Carlo  method  is  presented  as  an  approximation  of 
the  Neumann  series  solution  of  the  Integral  transport  equation. 

A  discussion  of  the  method  developed  to  bias  the  sampling  from  the 
source  angular  distributions  is  given,  and  a  description  of  the 
Russian  Roulette  and  splitting  techniques  utilized  by  the  procedure 
is  presented. 
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A  comparison  of  results  obtained  with  the  procedure  and  from 
experimental  data  taken  in  the  duct  systemizatlon  study  is  shown 
to  confirm  the  validity  of  the  procedure. 

The  ability  of  the  procedure  to  treat  multiple  scattering 
rigorously  and  the  versatility  offered  by  the  procedure  in  geometric 
description  allow  a  fairly  accurate  analysis  of  neutron  and  gamma- 
ray  radiation  streaming  through  a  multibend  duct.  Use  of  the  pro¬ 
cedure  should  lead  to  a  more  confident  prediction  of  the  flux 
streaming  through  multibend  ducts  than  has  been  obtainable  with 
some  of  the  less  exacting  methods. 
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I .  INTRODUCTION 


The  methods  presently  available  for  calculating  the  amount  of 
radiation  that  penetrates  a  multibend  duct  lack  the  versatility 
that  Is  necessary  to  give  a  complete  analysis  of  the  multiple 
scattering  which  Is  certain.  In  many  geometric  configurations,  to 
be  responsible  for  a  significant  portion  of  the  flux  penetrating  the 
duct.  Two  of  the  methods  most  commonly  used  -  the  single-scattering 
method  and  the  albedo  method  -  are  described  In  Reference  1,  These 
methods  have  given  good  results  for  a  certain  class  of  duct  config¬ 
urations,  namely,  those  that  have  a  high  length-to-dlameter  ratio. 

As  pointed  out  In  Reference  1,  however,  these  methods  have  failed 
to  predict  experimental  data  for  ducts  with  low  length-to-dlameter 
ratios.  In  References  2  and  3»  applications  of  diffusion  theory  to 
the  problem  of  radiation  penetrating  shields  containing  multibend 
ducts  are  described,  but  the  developments  of  the  theory  are  not 
supported  by  experimental  data. 

Because  of  the  success  of  a  recent  application  of  the  Monte 
Carlo  method  to  problems  Involving  the  scattering  of  neutrons  and 
their  penetration  through  straight  cylindrical  ducts  (Ref.  4),  It 
was  decided  to  modify  the  procedure  so  that  raultlbend  ducts  might 
be  considered.  Routines  providing  for  the  treatment  of  gamma-ray 
scattering  and  for  the  splitting  and  Russian  Roulette  techniques 
were  also  added  to  the  procedure. 
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This  modified  version,  referred  to  as  The  Multlbend-Duct  Pro¬ 


cedure  (LO5),  has  been  coded  for  the  IBM  7090  for  use  in  the 
analysis  of  neutron  or  gamma-ray  radiation  streaming  through  a 
m.ulfcibend  duct.  The  procedure  is  designed  to  calculate  the  direct- 
beam,  scattered  fluxes  and  doss  rates  at  from  1  to  30  detector 
positions  at  points  of  interest  within  or  near  a  multlbend-duct 
configuration.  Print  options  are  provided  so  that  the  energy  and 
angular  distributions  of  the  scattered  flux  at  each  detector  posi¬ 
tion  may  be  obtained. 

One  of  the  greatest  assets  of  this  procedure  is  that  multiple 
scattering  within  a  multlbend-duct  geometry  may  be  treated  rigor¬ 
ously.  The  procedure  allows  an  accurate  description  of  most 
multibend  ducts .  It  is  anticipated  that  the  multlbend-duct  pro¬ 
cedure  will  prove  to  be  a  valuable  tool  in  the  analysis  of  experi¬ 
mental  data  and  in  determining  the  validity  and  range  of  applica¬ 
bility  of  some  of  the  simpler  methods  developed  to  calculate  the 
radiation  penetrating  a  multibend  duct. 

A  description  of  the  geometry  and  of  the  Monte  Carlo  method 
utilized  by  the  procedure  is  given  in  Section  II.  An  evaluation 
of  the  procedure,  based  on  comparisons  of  calculated  and  mea¬ 
sured  data,  is  given  in  Section  III. 
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II. 


METHOD 


Tne  procedure  Involves  the  application  of  the  Monte  Carlo 
method  to  compute  the  energy  and  angular  distributions  and  the 
Intensity  of  the  scattered  neutron  or  gamma  flux  that  penetrates 
a  multibend  duct.  The  unscattered  flux  Is  also  computed  and 
recorded  separately  so  that  the  total  flux  at  any  detector  point, 
as  well  as  the  magnitudes  of  the  scattered  and  unacattered  f lux, ' 
may  be  determined.  Options  are  provided  for  splitting  and  Russian 
Roulette  and  for  biasing  the  source  angular  distributions. 

2. 1  Geometry  and  Source  Description 

The  geometrical  description  of  a  multlbend-duct  configuration 
is  accomplished  In  the  manner  described  In  this  subsection.  A  base 
coordinate  system  is  chosen  with  its  origin  lying  outside  of  the 
geometry  to  be  described  (Fig,  l).  The  multibend-duct  configuration 
is  then  divided  Into  sections  by  a  set  of  planes,  so  that  between 
any  two  consecutive  planes  the  duct  contains  no  bends  and  has  a 
constant  diameter.  Associated  with  the  second  of  any  two  consecu¬ 
tive  planes  Is  a  set  of  radii  consisting  of  the  outer  radii  of  each 
cylindrical  shell  region  between  the  two  planes.  In  tracking 
particles  as  they  scatter  through  this  geometry,  a  double  subscript 
la  used  to  Identify  the  regions;  for  Instance,  2p^  would  be  the 
macroscopic  cross  section  for  the  1th  cylindrical  region  associated 
with  the  pth  plane.  In  the  base  coordinate  system,  the  pth  plane 
Is  defined  by  the  point  of  Intersection  (Xp,Yp,Zp)  of  the  plane 
and  its  normal. 
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NPC  17,288 


Planes 


1.  Planes  are  identified  by  a  single  point 
(Xp,  Yp  ,  Zp)  in  the  base  coordinate  system 

2.  Coordinates  of  source  points, ,  and 
detector  points,  S  ,  must  be  given  in 
terms  of  the  base  coordinate  system 

3.  Beginning  v^tl^tjie  second  plane,  a  set  of 
3  vectors  (A,  B,C)  is  input  to  describe  the 
transformation  from  the  base  coordinate  system 
to  the  prime  coordinate  system  that  has  the  Z' 
axis  coincident  with  the  axis  of  the  duct  on 
the  origin  side  of  the  plane  (see  Fig.  2) 


Figure  1.  Mulfibend  Duct  Geometry 
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A  prime  coordinate  system  is  also  associated  with  each  plane 
(beginning  with  the  second).  The  Z'  axis  of  this  system  is  coin¬ 
cident  with  the  line  of  symmetry  of  the  regions  between  the  two 
consecutive  planes.  The  transformation  from  the  base  coordinate 
system  to  the  prime  coordinate  system  is  accomplished  by  three 
vectors:  A  =  a^,  a2,  the  vector  to  the  origin  of  the  prime 

coordinate  system;  B  r  b^^,  b2>  b^^  the  vector  to  the  point  (0,  0,  10) 
in  the  prime  coordinate  system;  and  C  =  c^^,  cg*  c^,  a  vector  to  the  - 
point  (10,  0,  O)  in  the  prime  coordinate  system.  Therefore,  the 
relationship  between  a  point  (X',  Y',  Z')  in  the  prime  coordinate 
system  and  the  same  point  (X,  Y,  Z)  in  the  base  system  is  as 
follows : 

Letting  X  s  X,  Y,  Z  be  a  vector  to  the  point  (X,  Y,  Z), 

then  X'  -  (X  -  A)  •  (C  -  A) 

10  " 


,  _  (X  -  A)  .  [(B  -  A)X  (C  -  A)] 


and 


Y'  = 


Z' 


100 

(X  -  A)  ♦  (B  -  A) 
10 


The  angles  k'  and  0'^  which  define  the  particles  direction 
at  the  point  (X,  Y,  Z)  or  (X' ,  Y',  Z')^are  given  by 


k'  -  cos 


-1 


(B 

>  1 

M  1 

1 

B  -  A  1 

cos  (if  =  <°  -  *)  •  i 


«  a. 

C  -  A 

Sin  k' 

and 
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Figpra  2.  Transformation  of  Coordinate  System 


sin  r  .  I(g  -  A)  X  (c  -  A)]  ■  E _  _ 

|(B  -  a)  X  (C  -  A)|sln  k' 

where  I  Is  a  unit  vector  In  the  direction  of  the  particle. 

A  set  of  from  1  to  30  source  points  may  be  used  to  approxi¬ 
mate  a  neutron  or  gamma-ray  source.  The  source  energy  may  be 
monoenergetlc  or  an  arbitrary  spectrum,  but  the  spectrum  must  be 
the  same  for  all  source  points.  There  may  be  a  different  angular 
distribution  for  each  source  point,  and  a  different  importance 
function  may  be  used  to  sample  from  each  of  the  angular  distribu¬ 
tions  . 

2.2  Scattering  Theory 

The  Monte  Carlo  method  has,  in  general,  been  considered  as 
a  method  of  solving  Integral  equations.  The  method  lends  Itself 
well  to  those  multiple  integrals  that  have  their  limits  complicated 
by  the  boundary  conditions  imposed.  Thus,  it  has  been  used  quite  ex¬ 
tensively  in  the  solution  of  neutron  and  gamma-ray  transport  prob¬ 
lems  involving  complex  geometries. 

The  method,  as  presented  here,  approximates  the  Neumann  series 
solution  of  the  Integral  transport  equation: 

-  _ 

P(R)  =  Y  In  (R), 

n=o 

where  F(R)  represents  the  total  flux  at  detector  position  R,  and 
In(R)  is  the  contribution  to  the  flux  from  the  order  of  scattering. 
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If  the  series  P(r)  converges,  it  may  be  approximated  by 

_  ^max  _ 

F(R)  =  Z  In(R),  (1) 

n»0 

where  the  sum  of  the  contributions  to  the  flux  from  orders  of 
scattering  above  nf^^^  negligible  when  compared  to  the  sum  of 
the  contributions  from  orders  of  scattering  less  than  or  equal  to 

^max» 


The  direct  beam,  or  unscattered  flux,  at  position  R  Is  given 
by  the  first  term  of  the  series, 

•  (R-Ro) 

J  Zrp(Ro-t  sn',Eo)  ds 


Ir(R)  »  - ; - To - 


- 2 

R“Rn 


6  (n<-no)f 


where  S(Roj.nosEo)  Is  the  source  strength  at  the  position  Rg 

expressed  in  the  terms  of  neutrons  per  second 
per  steradlan  In  the  direction  Hq  with  energy 

Eoj 

sCfi'-nQ)  Is  the  Dirac  delta  function,  and 

n*  Is  a  unit  vector  In  the  direction  R-Rg* 

The  Integral  expression 

•  (R-Ro) 

J'  2t(Ro  +  ^ j  ^o  ) 

0 


represents  the  number  of  mean-f ree-path  lengths  between  Rg  and  R 
and  Is  written  in  this  form  to  indicate  that  the  macroscopic  cross 
section,  Zt,  is  a  function  of  position. 
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The  first  term  of  the  series,  Io(R),  may  be  determined  exactly, 
but  the  rest  of  the  terms  are  approximated  by  the  Monte  Carlo  method. 
A  monoenergetic  point  source  will  be  assumed  for  the  presentation 


of  the  method.  The  single-scattered  flux  is  then  represented  by 


Ii(R) 


Zt(Ri,Eo)  exp  -  f  Zt(Ro+^o^Eo)  ds 

El-^o  *  L  o 


-  _  r  _  dRi 

W(Ri,Eo)  f(no,Eo - ►il',E')  exp  -  J  Zt(Ri  +  '  ;E' )  ds 


•  (R-Ri) 


where  ZtCRi^Eq)  Is  the  macroscopic  total  cross  section  at  position 
R]_  and  at  energy  Eq, 

is  the  ratio  of  the  macroscopic  scattering 
to  the  macroscopic  total  cross  section  at 
position  R^  and  energy  Eq, 

f(i2Q,Eo — ^fl',E')  is  the  probability  density  function*  for  scattering 

from  direction  to  I^n  > 

dRi  =  Ri^  sin  kg  dkgdjZfgdRi, 

where  Rj  “| Ri-Rq | » 

kg  is  the  polar  angle  In  the  direction  Q,q,  and 
0Q  is  the  azimuthal  angle  in  the  direction  Hq. 


*A  detailed  discussion  of  this  function  along  with  methods  of 
sampling  scattering  angles  from  the  function  is  given  in 
Section  2.3  of  this  report. 
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In  general,  for  n  2:  1, 


I 


n 


S(Ro,fto>Eo)  K(Rn_i,Rn)  W(Rn,En_i) 


En(R,n',E')  dRn .  dR^, 

where  _  _  _ 

n-1 


0 


K(Rn-i,Rn)  =  IT  2t(Ri+i 

1=0 


E^ )  exp 


2ip(Rj^"f'' Ej[ )  ds 


^  (^1  _  Ej[_ ,  Ej_ )  _  __ 

%+!"% 


(2) 


(3) 


with 


- ►no,Eo)  =  1  for  1  =  0, 


W(Rn,En-i) 


n 


1=1 


2s(Ri>Ej^_2) 

Et(^1> Ei_l) 


(M 
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and 


E„(R,  E')  =  exp 


«'  •  (R  -  Rn) 

r 2T(Rn  +  S  a  E' )  ds 


1 

«n-l»  En.i - -ii',  E'  )  I  -  _  2  *  (5) 

The  Monte  Carlo  method  utilized  In  approximating  Equation  1 
is  a  history-generating  process  whereby  each  history  gives  an 
estimate  of  each  term  of  the  series  but  the  first,  and  thus  gen¬ 
erates  an  estimate  of  the  scattered  flux.  One  history  is  the 
process  of  sampling  path  lengths  and  scattering  angles  from  the 
expressions  K(Rjj_2^,  R^^)  and  evaluating  the  estimators 
E|n(R,  E')  and  W  (R^,  ^n-1^  first  term  of  the 

series.  It  will  be  noticed  that 


K(R„.i,  R„)  =  K(R„.2,  R„, 


)  S 


■'«n>  ^n-l) 


exp 


^  n-l  ‘  (^n  ■  ^n-1^ 

J  ^  *  s  a  „.i,  E„.j)  ds 


f(  fl 


n-2*  n-2 


•“n-l'  ' 


SO  that  the  position  of  the  nth  collision  is  dependent  upon  the 
position  of  the  n-lst  collision.  Thus,  samples  drawn  from 
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K  (Rn-2'  ^n-1^  obtain  an  estimate  of  Ij^_2(R)  may  be  reused  In 
obtaining  an  estimate  of  In(R). 

Evaluating  the  expressions  Ej^(R,  0,  E' )  and  W  ^h-l^ 

for  the  values  of  R^  sampled  from  KCR^.i*  R^)  Slves  an  estimate 
of  the  magnitude  of  the  flux  at  R,  which  has  direction  '  and 
energy  E' .  Therefore,  by  allocating  storage  locations  for  energy 
and  angle  groups,  an  estimate  may  be  added  Into  the  proper  groups; 
and,  after  a  sufficient  number  of  histories  are  processed,  histo¬ 
grams  of  the  energy  and  angular  distributions  of  the  scattered 
flux  can  be  determined.  Estimates  of  the  total-scattered  flux 
at  R,  S^(R)  are  obtained  by  adding  the  estimates  from  each  order 
of  scattering. 

The  final  estimate  of  the  scattered  flux  at  a  position  R  Is 
taken  to  be  the  average  value  of  the  estimates  obtained  from  the 
Individual  histories : 


S(R) 


H 


I 


J 


where  S(R)  Is  the  scattered  flux  at  the  position  R,  as  estimated 
by  the  Monte  Carlo  method,  and  Sj^(R)  represents  the  estimates  of 
the  scattered  flux  obtained  from  the  individual  histories.  An 
estimate  of  the  variance  of  the  fluxes  from  the  Individual 
histories  is  given  by 


H 

Z 

irl 


2 


W 
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When  a  sufficient  number  of  histories  have  been  processed. 


'^max 

S(R)  ^  2  In(R) 

n“l 

and 

IqCr)  +  S{R)  =  F(R). 

2. 3  Sampling  Methods 

In  the  Monte  Carlo  method  utilized  by  this  procedure,  path 
lengths  and  scattering  angles  are  chosen  by  a  random  process  to 
obtain  the  distribution  of  the  collision  points  in  a  geometric 
configuration  containing  a  multibend  duct.  The  concept  of  weight¬ 
ing  is  also  Introduced  to  allow  sampling  from  other  than  the  actual 
physical  distributions  of  source  angles  and  collision  densities, 
and  to  adjust  for  the  fact  that  no  absorption  collisions  are 
allowed. 

A  method  of  biasing  the  sampling  from  the  angular  distribu¬ 
tion  of  a  polnt-lsotroplc  or  a  point-anisotropic  source  was  deve¬ 
loped  so  that  those  angles  which  give  the  source  particle  a  higher 
probability  of  reaching  the  detector  may  be  given  more  emphasis  than 
they  would  ordinarily  receive  in  straightforward  sampling.  For 
this  biased  sampling  scheme,  the  unit  sphere  about  the  source 
point  is  divided  into  areas  which  may  be  unequal  in  size  but  are 
considered  to  be  of  equal  importance.  The  divisions  in  area  are 
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made  with  respect  to  the  polar  angle,  giving  I  angle  sectors  that 
are  bounded  by  I  -  1  polar  angles.  Each  sector  is  considered  to 
be  of  equal  Importance,  so  that  the  total  number  of  histories  to  be 
run  are  divided  equally  among  the  sectors.  The  adjustment  to  the 
weighting  factor  for  the  1^  vector  to  remove  the  bias  introduced 
by  the  sampling  process  Is 


W 


1 


I 


n(©)  sin  9d9 


n(9)  sin  9d0 


(6) 


where  I  is  the  total  number  of  sectors,  and  are  the  polar 

angles  bounding  the  1^  sector,  and  n(9)  Is  the  angular  density 
function  describing  the  angular  distribution  of  particles  emitted 

from  the  source , 

The  cosines  of  the  polar  angles  defining  the  sectors  are 
listed  as  input  for  the  procedure,  and  the  weighting  factors  of 
each  sector  are  listed  when  the  angular  distribution  is  anisotropic. 
For  the  case  of  an  isotropic  distribution,  the  weight  adjustment 
factors  are  calculated  by  the  machine  by  use  of  Equation  ' S, ,  which 
reduces  to 


W.  =  I 


cos  -  cos  kj 


when  n(9)  is  the  isotropic  distribution. 
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The  patn  lengths,  ®|ki  “  a-i’e  sampled  from  the 

distribution  function 


dRi, 


1^1 

r 

RN  =  I  2<;p(RjL?E.i_]^)  exp 

-  J 

0 

where  RN  is  a  random  number,  and 

%-l  '  (Rl-Ri-i) 


ZT(Ri,Ei_l)  exp 


2(p{Rj^_l  ^l-l>®l-l) 


is  one  of  the  integral  expressions  given  in  Equation  3. 

Pi  =  Pi-l-h  P-1^1-1  gives  the  position  of  the  1th  collision.  After 
determining  the  position  of  each  collision,  the  weighting  function 
is  multiplied  by  Ss(Riji )/Et(Ri^ ),  since  only  scattering 
collisions  are  allowed.  Prom  each  collision  point  an  estimate  of 
the  flux  reaching  each  of  the  detector  points  is  made,  assuming 
uhat  no  further  collisions  occur  between  that  point  and  the 
detectors.  The  estimating  function  is 

W(Rn.En_i)  En(R,ft',E' ), 

as  given  in  Equations  4  and  5. 

A  close  examination  of  the  expression 

exp 

o 


ft'  •  (R-Rn) 

/  ZT(Rn  +  Sft',E'  )  ds 
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contained  within  En(R,^^ ' /E*  ),  reveals  that  collisions  occurring 
near  a  detector  point  will  contr'ibute  much  more  to  the  scattered 
riux  than  those  occurring  farther  away.  This  being  the  case,  it 
is  advisable  to  ensure  that  a  sufficient  number  of  collisions  are 
sampled  in  those  areas  near  the  detector.  In  order  to  improve  the 
sampling  of  collisions  in  areas  near  a  detector,  an  option  pro¬ 
viding  for  the  splitting  and  Russian  Roulette  techniques  is  in¬ 
cluded  within  the  procedure.  Splitting  and  Russian  Roulette  re¬ 
sult  in  a  change  in  the  number  of  particles  being  followed.  Con¬ 
sequently,  the  particle  weights  must  be  correctly  adjusted  to 
eliminate  the  bias  introduced  by  these  techniques.  The  option  re¬ 
quires  that  a  splitting  factor,  M,  and  a  set  of  splitting  bound¬ 
aries,  Hg,  be  listed  as  input.  The  factor,  M,  is  an  Integer  and 
determines  the  number  of  particles  into  which  a  single  particle 
will  be  split  or  the  number  of  particles  that  will  be  combined 
into  one  particle  upon  crossing  a  splitting  boundary.  The  split¬ 
ting  boundaries  H2,  ...  ...  are  listed  in  ascending 

order  and  define  the  loci  of  points  lying  H^,  E2t  Hg,  ..., 

relaxation  lengths  from  the  detector  point.  The  quantity 

ii'  •  (R  -  Rj^) 

ll  -  j  Z^(R„  f  Sii  ',  E')ds 
o 

in  Equation  5  gives  the  number  of  relaxation  lengths  between  the 
nth  collision  point  and  the  detector;  therefore  B,  the  index  of  the 
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largest  value  of  Hg<H,  Is  the  number  of  splitting  boundaries  be¬ 
tween  the  collision  point  and  the  detector.  If  B'  and  B  are 
the  number  of  splitting  boundaries  between  the  detector  point  and 
the  n-lst  and  the  n;^  collision  point,  respectively,  then  |B'-Bj 
is  the  number  of  splitting  boundaries  that  a  particle  crosses  In 
going  from  the  n-ls_t  to  the  n^  collision  point.  If  B'>B,  the 
particle  Is  moving  towai^d  the  detector  and  splitting  occurs.  In 
this  case,  the  particle  Is  split  Into  particles,  so  that  the 

weight  of  each  particle  becomes 


W 


n-1 


_  w' 


n-1 


M- 


,B-B' 


where  was  the  weight  of  the  original  particle  before  split¬ 

ting.  If  B'<B,  the  particle  is  moving  away  from  the  detector  and 
the  Russian  Roulette  technique  is  applied.  A  random  number,  RN, 

Is  generated,  and  If  RN>(j^)  ,  the  tracking  of  the  particle 

Is  terminated;  but  if  RN<(^)®“®',  the  weight  Is  multiplied 

I  ^ 

by  (M)  and  tracking  of  the  particle  Is  continued.  A  proper 
selection  of  the  splitting  factor  and  a  proper  placement  of  the 
splitting  boundaries  should  improve  the  sampling  of  the  collision 
points  in  the  areas  near  the  detector. 

After  an  estimate  of  the  flux  reaching  each  detector  from  the 
nth  collision  has  been  made,  the  procedure  returns  to  the  loca¬ 
tion  of  the  nth  collision  to  choose  a  scattering  angle, 
ij/  s  cos  (  *  ^n-1^* 
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For  gamma  rays,  the  scattering  angle  Is  chosen  from  the 
Kleln-Nishina  differential  scattering  cross  section  for  Compton 
scattering  or  from  an  isotropic  distribution  for  the  pair-production 
process.  The  pair-production  process  is  assumed  to  be  a  scatter¬ 
ing  process  in  which  two  0.50-Mev  photons  are  given  off  in  opposite 
directions.  The  procedure  chooses  the  scattering  process  by  com¬ 
paring  a  random  number  with  the  ratio  of  the  Compton-scatterlng 
cross  section  to  the  sum  of  the  Compton-scatterlng  and  pair- 
production  cross  sections. 


RN 


pe  ^c^^n-1^ 

Pe(f^n)  M-c^^n-l^  “  ^  pp^*^*  ^n-1^ 

J 


where  Pg(Rj^)  is  the  electron  density  at  the  position  R^  , 
is  the  Compton-scatterlng  cross  section,  p  j(Rn^ 
is  the  atomic  density  of  the  element  j  at  the  position  Rj^,  and 
)ipp(j,  Efi-l^  pair-production  cross  section  for  element  J. 

The  summation  is  over  all  elements  J  at  the  position  If  RN 

is  less  than  or  equal  to  the  right-hand  member  of  Equation  6, 
the  scattering  process  is  taken  to  be  Compton  scattering.  If 
RN  is  greater  than  the  right-hand  member  of  Equation  6,  the 
scattering  process  is  taken  to  be  the  pair-production  process. 

The  angular  distribution  of  the  photons  produced  by  the  pair- 
production  pivDcess  is  assumed  to  be  isotropic  in  the  laboratory 
system,  so  that  ^ ,  the  polar  angle  of  emission  with  respect  to 
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the  previous  direction  of  flight  Q.  Is  given  by  the  following 

equation : 

^  =  cos-1  (2RN-1), 

where  RN  Is  a  random  number.  To  prevent  the  tracing  of  Individual 
histories  for  both  of  the  photons  given  by  the  pair- production 
process,  both  photons  are  assumed  to  be  emitted  In  the  same  direc¬ 
tion  and  are  combined  Into  one  photon  having  twice  the  weight  of 
the  original  photon.  The  bias  Introduced  by  assuming  that  both 
photons  are  emitted  In  the  same  direction  will  automatically  be 
eliminated  by  running  a  large  number  of  histories. 

Sampling  of  the  scattering  angle  from  the  Kleln-Nlshlna 
differential  scattering  cross  section  for  Compton-scattered  pho¬ 
tons  Is  done  by  the  rejection  technique  described  on  the  bottom 
of  page  65  In  Reference  5» 

For  neutrons,  the  sampling  of  the  scattering  angle  Involves, 
first,  a  choice  of  the  target  element  and,  then,  a  choice  between 
the  elastic- and  Inelastic-scattering  processes  for  that  element. 
The  target  element  J  Is  selected  so  that 

^n-l)  ^  'T  ^n-l^ 

..  _ - _ - _  <  RN  <  ^  _  > 

k=l  E„_i)  k=l  2t(R„,  E„.i) 

where  P  Is  the  atomic  density  of  the  element  k  at 
aij(k,  Is  the  microscopic  total  cross  section  for  the 
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element  k,  and  2  macroscopic  total  cross  sec¬ 

tion  at  the  position  R^.  The  differential  elastic-scattering 
cx‘oss  section  or  the  differential  Inelastic-scattering  cross  sec¬ 
tion  for  the  element  J  is  taken  to  be  the  distribution  of  the 
scattering  angle  f  ,  depending  upon  whether  a  random  number  RN  Is 
less  than  or  greater  than 


a 


a 


el 


el 


(J,  E 


n-1 


(Jj>  Ej^_]_ 


) 

n 


bin(Jji  ®n-l^ 


The  procedure  assumes  that  the  inelastic  scattering  process  Is 
isotropic  In  the  center-of-mass  system  for  all  elements.  The 
center-of-mass  scattering  angle,  Is  then  determined  by  the 

equation 


-  cos"^  [  2RN  -  l] 

and  f  ,  the  angle  In  the  laboratory  system  corresponding  to  ^ 
Is  determined  by 

_l  1  A'  cos  ^  ' 

^  z  cos - 1 - -  , 

_  (1  +  2  A'  cos  ^  '  +  A'^)2_ 


where 


A' 


1 


(Aj  +  l) ^ 


X 

2 
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with  Aj  being  the  atomic  weight  of  the  element  J  and  with  £  being 
the  excitation  level  of  the  target  nucleus.  The  excitation  levels 
of  the  target  elements  are  selected  from  a  table  of  input  data 
listing  the  probabilities  of  exciting  the  various  levels  of  the 
target  nucleus  as  a  function  of  the  Incident  neutron  energies. 

Routines  are  Included  In  the  procedure  for  the  choice  of  ^ 

In  the  special  cases  when  the  differential  elastic-scattering 
cross  section  may  be  considered  to  be  Isotropic  in  the  laboratory 
system  or  isotropic  In  the  center-of-mass  system.  The  procedure 
also  allows  for  a  table  of  Input  values  of  the  quantity 

^  sin  V'’  dK/’ 

IT 

f  ^  sin  f  '  d  ^  ' 
j  dfi 
o 

for  elastic  scattering  as  a  function  of  energy,  angle,  and  element. 
If  these  tables  are  used,  the  set  of  values  for  the  energy  In  the 
.table  nearest  the  energy  of  the  Incident  neutron  Is  chosen,  and  a 
linear  interpolation  Is  performed  to  determine  the  value  of  ^ 
corresponding  to  a  random  number. 

Since  any  azimuthal  angle  ^  about  the  previous  direction 
of  a  particle  Is  equally  probable,  random  values  of  (|)  may 
be  found  by  evaluating  the  Integral  in  the  following  equation  for 
C|)  ,  where  RN  Is  a  random  variable: 


0 


/ 
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The  directional  vectors  may  be  defined  with  two  components 
k^,  the  polar  angle  with  respect  to  the  Z  axis,  and  0,  the  angle 
between  the  projection  of  the  vector  In  the  (x,y)  plane  and  the  x 
axis.  The  components  kn  and  0n  of  the  vector  determined 

from  the  scattering  angles  k^_i,  and  j2Cn_i  by  the  following 

spherical  triangle  relationships: 


kn  = 
cos  ^n  “ 
sin  ^n  = 


cos"^j^cos  kn_i  cosf  +  sin  kn_i  sln^  cos 
cos  cos  A0^  -  3lnj2^n-l  ^^oA0n>  and 

sin  0YI-1  cosA^n  cos^n-l  sinA0j^> 


where 


and 


cos  i/  -  cos  cos 

cosAjZ^n  =  sin  kn_i  sin  kn 


sinAp'ri  = 


sin  Tj!  sin 
sin  kn 
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III.  COMPARISON  WITH  EXPERIMENTAL  DATA 


General  Dynamic  s/Fort  VJorth  is  presently  engaged  in  a  duct 
systemizatlon  and  penetration  study  (Ref,  6).  Experimental  data 
for  several  straight  cylindrical  ducts  of  different  lengths  and 
diameters  have  been  obtained  in  support  of  the  study^  and  it  is 
planned  to  continue  the  experiments  with  multlbend-duct  configura¬ 
tions.  A  literature  survey  for  other  experimental  data  Involving 
multibend  ducts  revealed  that,  of  the  data  available,  none  was 
reported  in  sufficient  detail  to  permit  conclusive  comparisons 
with  Monte  Carlo  results. 

Since  the  Multibend- Duct  Procedure  may  also  be  applied  to 
straight  cylindrical  ducts  and  since  the  only  difference  between 
the  handling  of  straight  and  multibend  ducts  by  the  procedure  is 
the  difference  in  the  transformations  used,  it  was  decided  to  use 
the  experimental  data  taken  with  straight  cylindrical  ducts  to 
check  the  procedure. 

However,  multibend-duct  configurations  were  run  with  the 
procedure  and  checked  by  hand  calculations  to  ensure  that  the 
transformations  for  multibend  ducts  would  be  handled  properly. 

The  experimental  data  for  the  straight  cylindrical  ducts 
were  taken  by  inserting  aluminum  ducts  of  different  lengths  and 
diameters  through  a  porthole  in  the  bottom  of  a  6-foot  cubical 
tank  and  filling  the  tank  with  water  to  a  level  even  with  the  top 
of  the  duct  (Pig.  3).  Then  a  source  -  polonium-beryllium  for 
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Figure  3.  Experimental  Setup  for  Duct  Systemisatien  Study 
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neutrons  or  cobalt-60  for  gamma  rays  -  was  centered  at  the  mouth 
of  the  duct,  level  with  the  bottom  of  the  tank.  The  detector  -  a 
fast-neutron  dosimeter  (PND)  for  neutrons  or  an  anthracene 
scintillation  dosimeter  (ASD)  for  gammas  -  was  traversed  in  planes 
3  Inches  and  15  Inches  above  the  water  level.  Measurements  of  the 
dose  rates  for  the  3-inch  traverse  of  3-inch- and  6-inch-diameter, 
12-inch-long  ducts  were  chosen  for  comparison  with  results 
obtained  from  the  multibend-duct  procedure.  Figures  4  and  5 
show  the  comparison  of  the  gamma-ray  experimental  and  calculated 
data  for  the  3-inch  traverse  of  a  3-lnch-dlameter,  12-lnch-long 
duct  and  a  6-lnch-dlameter,  24-lnch-long  duct.  A  444-millicurie 
cobalt-60  source  was  used  in  the  experiment,  and  this  source  was 
assumed  to  be  Isotropic  in  the  calculations. 

The  comparison  with  experimental  data  is  somewhat  better 
than  was  expected.  Table  I  shows  the  statistical  variations  in 
the  calculated  scattered  dose  rates,  and  they  became  quite  large 
as  the  detector  point  was  moved  farther  from  the  centerline. 

There  was  also  some  question  as  to  how  well  the  ASD  response 
curve  and  the  flux-to-dose  conversion  curve  used  in  the  calcu¬ 
lations  would  agree  over  the  energy  range  considered.  The  energy 
range  Itself  was  in  questlon_,  since  in  the  calculations  all  photons 
with  energies  below  0.2  Mev  were  assumed  to  be  absorbed.  Later  it 
was  found  that  this  was  a  questionable  assumption,  because,  even 
though  the  ASD  response  to  photons  below  0.2  Mev  is  fairly  low, 
there  is  a  good  possibility  that  if  the  spectrum  of  the  scattered- 


32 


Gamma-Ray  Dose  Rate 


NPC  17,292 


Distance  from  Centerline  (inches) 

Figure  5.  Comporison  of  Colculofod  and  Moaturod  Oammo*Roy 
Dos#  Ratot  for  o  6-(nch-Diomotor,  Rd-Inch-Long  Duct 
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photon  energies  at  the  detector  were  known  it  would  show  a  high 
peak  in  the  range  from  0.05  to  0.1  Mev.  Measurements  of  the  gamma- 
ray  spectra  from  cobalt-60  gamma  rays  scattered  by  1-  to  12-inch- 
thick  polyethylene  slabs  show  a  high  peak  at  0.1  Mev  for  the  l-inch 
slab,  which  gradually  shifts  to  O.05  Mev  as  the  slab  thickness 
increases  to  12  inches  (Ref.  7)«  Polyethylene  and  water  have  very 
similar  gamma- ray  attenuation  properties;  therefore,  peaks  in  the 
0.05-  to  O.l-Mev  energy  range  are  anticipated  in  the  scattered 
gamma-ray  spectra  at  the  ends  of  the  ducts  because  of  the  scattering 
of  the  cobalt-60  gamma  rays  within  the  water  surrounding  the  ducts. 
Probably  some  of  these  effects  tended  to  cancel  each  other  and, 
hence,  the  good  agreement  given  in  Figures  4  and  5  was  obtained. 

Figures  6  and  7  show  the  comparison  of  the  neutron  experi¬ 
mental  and  calculated  data  for  the  3-lnch  traverse  of  both  a  3-inch- 
and  a  6-inch-diameter  duet  12  inches  long.  A  comparison  with  data 
obtained  with  the  General  Electric  Flexible  Monte  Carlo  procedure 
FMC-N  (Ref.  8)  is  also  presented  in  Figure  7.  A  polonium-beryllium 
source  having  a  source  strength  of  3.98  x  10^  neutrons/sec  was  used 
to  obtain  the  experimental  data  for  neutrons.  The  source  angular 
distribution  was  slightly  anisotropic  because  of  the  cyllndrically 
shaped  capsule  containing  the  polonium-berylliiam.  The  anlsotrophy 
of  the  source  was  taken  into  account  in  the  calculation  of  the 
direct-beam  flux  at  the  detector  locations,  but  an  isotropic  source 
was  assumed  in  the  calculation  of  the  scattered  radiation. 

Table  II  shows  the  scattered  and  unscattered  components  of 
the  dose  rates  as  calculated  by  the  multibend  procedure. 
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Neutron  Dose  Rate  (mrem/Kr) 


NPC  17.293 


103 


Figure  6.  Comporison  of  Colculotod  and  Moaturod  Faff*Noufren 
Dose  Rotof  for  o  3-lnch-Diamotor,  IR-Inch^Long  Duct 
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The  measured  PND  response  curve  (Ref.  9)  shown  In  Figure  8  was 
used  in  the  calculations  rather  than  a  f lux-to-dose  conversion  factor. 
Therefore^,  any  dlsargreement  with  the  neutron  experimental  data  Is 
probably  due  to  a  statistical  variation  In  both  the  experimental  and 
calculated  data. 
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Neutron  Response 


NPC  17,295 
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Figure  8.  Fast-Neutron  Dosimeter  Response  Curve 
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IV.  CONCLUSIONS 


In  light  of  the  good  agreement  obtained  with  experimental 
data  for  both  the  neutron  and  gamma-ray  penetrations  through 
straight  cylindrical  ducts  and  under  the  assumption  that  the 
application  of  the  procedure  to  multibend  ducts  will  not  affect 
the  mechanics  of  the  calculations.  It  Is  concluded  that  the 
Multibend  Duct  Procedure  will  give  an  accurate  prediction  of 
neutron  or  gamma-ray  flux  penetrating  a  multibend  duct.  The 
ability  of  the  procedure  to  treat  multiple  scattering  rigorously 
and  the  versatility  offered  by  the  procedure  In  geometric  de¬ 
scription  allow  a  fairly  accurate  analysis  of  neutron  and  gamma- 
ray  radiation  streaming  through  a  multibend  duct.  This  should 
lead  to  a  more  confident  prediction  of  the  flvix  streaming  through 
suv.n  configurations  than  has  been  obtained  with  some  of  the  less 
exacting  methods. 

It  Is  anticipated  that  this  procedure  will  be  a  valuable 
tool  In  the  evaluation  of  the  validity  and  the  range  of  appli¬ 
cability  of  some  of  the  simpler  methods  developed  to  calculate 
the  amount  of  radiation  penetrating  a  multibend  duct. 
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